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Blazar candidates beyond redshift 4 observed with GROND 
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ABSTRACT 

The search for extremely massive high redshift blazars is essential to put strong constraints 
on the supermassive black holes formation theories. Up to now, the few blazars known to 
have a redshift larger than 4 have been discovered serendipitously. We try a more systematic 
approach. Assuming radio-loudness as a proxy for the jet orientation, we select a sample of 
extremely radio-loud quasars. We measure their black hole masses with a method based on 
fitting the thermal emission from the accretion disc. We achieve a precision of a factor of two 
for our measures, thanks to the observations performed with the Gamma-Ray Burst Optical 
Near-Infrared Detector (GROND). The infrared to optical GROND data allow us to observe 
directly the peak of the disc emission, thus constraining the overall disc luminosity. We obtain 
a small range of masses, that peaks at 10 9 ' 3 M Q . If some of our candidates will be confirmed 
as blazars, these results would introduce interesting constraints on the mass function of ex- 
tremely massive black holes at very high redshift. Moreover, all our blazar candidates have 
high accretion rates. This result, along with the high masses, opens an interesting view on the 
need of a fast growth of the heaviest black holes at very high redshift. 



Key words: Galaxies: 
frared: general 



active — quasars: general — radiation mechanism: thermal — in- 



1 INTRODUCTION 

Blazars are radio-loud Active Galactic Nuclei (AGN) with rela- 
tivistic jets directed at a small angle from our line-of-sight. It is 
convenient to define blazars those sources whose jet makes an an- 
gle # v to our line of sight comparable to or smaller than 1/T, 
where V is the bulk Lorentz factor of the emitting region of the 
jet. At these angles the jet emission is strongly boosted by rela- 
tivistic effects. Hence, blazars should be in principle visible even 
at very high redshift. Nevertheless, up to now, only few of them 
are known at z > 4, and they have all been discovered serendip- 
itously. The typical blazar Spectral Energy Distribution (SED) is 
dominated by the non-thermal jet emission, and is characterized by 
two non-thermal humps, respectively attributed to synchrotron (at 
low frequency) and inverse Compton (high frequency) processes. 
Generally, the emission from the other structures of the blazar is 
not or barely visible, because of this jet dominance (i.e. the thermal 
emission from the accretion disc, that one should expect to observe 
in the optical-UV band, is totally or partially covered by the syn- 
chrotron component). As the bolometric luminosity increases, the 
two humps shift to lower frequencies (for a different interpreta- 
tion of this effect, see Giommi et al. 2012). Indeed, the synchrotron 
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component can shift enough to let the thermal emission from the 
accretion disc to emerge. In this case, the blazar would show to the 
observer an optical spectrum with features completely similar to 
those displayed by a radio-quiet AGN with an analogous accretion 
disc, not contaminated by other components. To illustrate this point, 
in Fig. [TJ we show the optical spectra of two quasars at the same 
redshift: one is radio-quiet, and the other is strongly radio-loud. 
The two spectra are indistinguishable. This is because at z > 4 
we can only see the most powerful objects, and there is a correla- 
tion between the luminosity of the accretion disc, the jet power and 
the peak frequencies of the non-thermal emission (Ghisellini et al. 
2011, Sbarrato et al. 2012a, Ghisellini et al. 2010a,b). For larger 
luminosities, the synchrotron peak shifts to sub-mm frequencies, 
leaving the accretion disc component "naked" and thus observable. 
An analogous effect occurs for the high-energy hump, whose peak 
is shifted below 100 MeV, making its detection by the Large Area 
Telescope (LAT) onboard the Fermi Gamma-Ray Space Telescope 
(Atwood et al., 2009) unlikely. Instead, a powerful blazar can be 
more easily detected at frequencies just before the inverse Comp- 
ton peak by hard X-ray instruments, e.g. the Burst Alert Telescope 
(BAT) onboard the Swift satellite (Gehrels et al. 2004). This effect 
is enhanced by the high redshift itself, therefore it is surely more ef- 
fective to use a hard X-ray telescope instead of a 7-ray instrument 
to detect high power 2 > 4 blazars. 
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Figure 1. Comparison of the optical spectra of an extremely radio-loud 
(upper panel) and a radio-quiet AGN (lower panel) with the same redshift, 
both from the sample of S 1 1 . Note that there are no differences between the 
two spectra that could suggest their difference in radio-loudness. 



By measuring the power of the source close to where it peaks, 
X-ray data are also crucial to classify the source as a blazar, since 
misaligned or even slightly misaligned jets correspond to much less 
luminous sources than their aligned counterparts. Following these 
guidelines, we could recently classify B2 1023+25, a z ~5.3 radio- 
loud source, as a blazar, thank to its hard and strong X-ray emission 
(Sbarrato et al. 2012b). 

The importance of observing blazars at z > 4 resides in the 
fact that we expect their black holes (BH) to be very massive, since 
we observe only the most powerful objects: A/bh can even exceed 
a billion solar masses. 

A large number of extremely massive BHs at z > 4 could 
potentially put strong constraints on the study of their formation. 
Because of their peculiar orientation, they have a great statistical 
relevance: for each observed blazar, in fact, there should be 2F 2 = 
450(r/15) 2 analogous radio-loud AGN with their jets directed in 
random directions (Ghisellini et al. 2010a; Volonteri et al. 2011) 
Discovering even a few blazars at high-z has a big impact on the 
study of early BHs, such as their comoving density at high-redshift. 
(Volonteri et al. 201 1; Rau et al. 2012). 

For these reasons, the search for high-z blazars can be as re- 
warding as the search for high-z radio-quiet sources, even if radio- 
quiet AGN are intrinsically more numerous than radio-loud ones. 
But to efficiently pursue this goal, we must search blazars in a sys- 
tematic way, instead of serendipitously as happened up to now. A 
systematic approach is possible if a complete quasar catalog is used 
as a starting point for searching for blazar candidates. In this work, 
we use the Quasar Catalog (Schneider et al. 2010) obtained from 
the 7 th Data Release (DR7) of the Sloan Digital Sky Survey (SDSS, 
York et al. 2000). 

Starting from this well defined (namely, flux limited and cov- 
ering a known portion of the sky) sample, we can select the best 
blazar candidates among the most radio-loud sources. In this way, 
once we will confirm the blazar nature of the candidates, we can 
construct the number density not only of the confirmed blazars, but 



also of their parent population (i.e. their misaligned counterparts). 
Moreover, by measuring the black hole mass in these systems, we 
can construct the number density of black holes of a given mass at 
high redshifts. 

In this work we will select a sample of blazar candidates cho- 
sen among the sample of SDSS quasars of Shen et al. (201 1, here- 
after Sll). The considered region of the sky is covered by the 
FIRST radio survey (1.4 GHz, with a flux limit of 1 mJy; White 
et al. 1997). We will select all quasars in that sample with z > 4 
and large radio-loudness. For each of them we will construct the 
IR-optical-UV spectral energy distribution (SED). All objects vis- 
ible from La Silla (Chile) have been observed with the Gamma- 
Ray Burst Optical Near-Infrared Detector (GROND; Greiner et al. 
2008), providing photometric data from the near IR to the optical. 
All objects have a SDSS spectrum, and all (but 4) are detected in the 
far IR by the Wide-field Infrared Survey Explorer (WISE; Wright 
et al. 2010). With the good wavelength coverage that is obtained, 
we are able to reliably estimate the black hole mass of the objects, 
by fitting the data with a standard accretion disc model. Along with 
the mass, we obtain the overall disc luminosity and, hence, the Ed- 
dington ratio for the whole sample. This work will provide a com- 
plete sample (i.e. flux limited) of possible blazar candidates, well 
described in their thermal emission and principal features, ^de- 
scribes the sample selection, shows the details of the GROND 
observations, ^explains the available data and the fitting method, 
in fj5]we discuss the results and their implications, while in Sj6]we 
summarize our work. 

We adopt a flat cosmology with Ho = 70 km s _1 Mpc -1 and 
Q M = 0.3. 



2 THE SAMPLE 

The SDSS provides the largest publicly available Quasar Catalog 
(Schneider et al. 2010), complete with both photometric and spec- 
troscopic data of 105,783 optically selected quasars. This catalog 
has been spectroscopically analyzed by Sll, who derived contin- 
uum and lines features, along with different virial mass estimates 
for all the objects. Sll cross-correlated the SDSS Quasar Catalog 
with the FIRST radio catalog, providing the 1.4 GHz radio fluxes 
and radio-loudnesses for all SDSS AGN with a FIRST radio coun- 
terpart. 

In order to select good high-redshift blazar candidates, we 
limit our research to the quasars with a radio matching and col- 
lect the most radio-loud quasars with z > 4. We obtain a sample 
of 3 1 quasars with z > 4 and radio-loudness 



R : 



fl.4GHz 



> 100. 



(1) 



The limit in radio-loudness should permit to preferentially select 
the objects with their jets oriented closer to our line-of-sight. This 
provides us a group of 31 possible candidates listed in Tab.[T]along 
with their most illustrative data. We include from Sll redshift, 
radio-loudness, radio flux, bolometric and CIV line luminosities 
and the best Black Hole mass estimate, where available. For those 
where the CIV line was visible, we include the overall Broad Line 
Region luminosity (£blr), derived following the template by Fran- 
cis et al. (1991). For the few objects already observed in the X-rays, 
we report in Tab.[T]the X-ray spectral index and flux (in the 0.3-10 
keV band). 

It is indeed worth to notice that among the 31 selected quasar, 
7 have already an X-ray detection. 2 of them can be classified as 
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blazars, thanks to their hard and intense X-ray fluxes (i.e. SDSS 
J083946.22+5 11202.8 and SDSS J15 1002.92+570243.3). Specif- 
ically, SDSS J15 1002.92+570243.3 is one of the firsts high red- 
shift (2 > 4) blazars ever discovered (Moran & Helfand, 1997). 
Moreover, we identified in a previous work our best candidate B2 
1023+25 as a blazar (Sbarrato et al. 2012b), thanks to the Swift X- 
ray data (Gehrels et al. 2004). 

Here, we want to estimate the central black hole mass (A/bh) 
of our blazar candidates. Since at 2 > 4 the virial mass method 
has to rely only on the CIV line, we use a different method to esti- 
mate A/bh, based on the direct fit of the accretion disc (see £ 14.21 ) . 
In order to fit it properly, we need a good coverage of the NIR- 
optical wavelength range. We used the spectroscopic and photomet- 
ric optical data provided by the SDSS DR7, along with the publicly 
available data from WISBj. To fill the gap between the wavelength 
ranges of these two releases, we observed our candidates from La 
Silla (Chile) with GROND, simultaneously in 3 NIR and 4 optical 
filters. The partial superposition of the GROND optical filters with 
the SDSS ones, allows to reduce the problems due to the possible 
variability of the objects. GROND data are very important, because 
we expect that the peak of the accretion disc emission falls between 
this wavelength range and the optical band provided by the SDSS. 
Since only 19 out of the 31 radio-loud high-2 quasars of the origi- 
nal sample are visible from La Silla, we reduce our sample to these 
19, listed in Tab.H These 19 sources have a very complete IR- 
optical coverage, and we focus on them to estimate the Mbh and 
also to compile a sample of good blazar candidates. 
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Table 2. Sources from the DR7 Quasar Catalog with z > 4 and radio- 
loudness R > 100, seen by GROND. Col. [1]: right ascension and declina- 
tion (i.e. SDSS name: "SDSS J. . . "); Col. [2]: redshift; Col. [3]: logarithm 
of the disc luminosity in ergs -1 , as derived in this work; Col. [4]: loga- 
rithm of the BH mass, as derived by SI 1 (solar masses); Col. [5]: logarithm 
of the BH mass as derived in this work (solar masses), with the confidence 
interval. 



3 GROND OBSERVATIONS AND DATA ANALYSIS 

The 7-band GROND imager is mounted at the 2.2m telescope of 
the Max-Planck-Gesellschaft (MPG), operated by MPG and ESO 
(European Southern Observatory) at La Silla (Chile). GROND is 
able to observe simultaneously in 7 filters, from the NIR K s (2400 
nm) to the g' band (450 nm), i.e. from 2200 to 360 nm. 

We carried out observations for all sources simultaneously in 
all 7 g',r',i',z',J, H, K s bands. The log of the GROND obser- 
vations and the related observing conditions are reported in Tab. 

El 

The GROND optical and NIR image reduction and photome- 
try were performed using standard IRAF tasks (Tody 1993), similar 
to the procedure described in Kriihler et al. (2008). A general model 
for the point-spread function (PSF) of each image was constructed 
using bright field stars, and it was then fitted to the point source. 
The fields of all sources in our sample are covered by the SDSS 
(Smith et al. 2002) survey, therefore the absolute calibration of the 
g',r',i' , z bands was obtained with respect to the magnitudes of 
SDSS stars within the blazar field. For all sources the J, H, K s 
bands calibrations were obtained with respect to magnitudes of the 
Two Micron All Sky Survey (2MASS) stars (Skrutskie et al. 2006). 

Tab.[3]reports the log of the GROND observations while Tab. 
[4]reports the observed AB magnitudes, not corrected for the Galac- 
tic extinction listed in the last column and taken from Schlafly & 
Finkbeiner(2011). 



Data retrieved from |http : / / irsa . ipac . caltech . edu/| 
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4 THE IR TO OPTICAL SED 

In the Appendix we show the IR to optical SED of all the 18 sources 
in our sample. We did not include SDSS J102623.61+254259.5, 
because we already applied to this source the accretion disc fitting 
method in our previous work (Sbarrato et al. 2012b). The first evi- 
dent feature of all the SEDs is the Lya forest, i.e. the prominent ab- 
sorption clearly visible in the SDSS spectra (black line) and in the 
g' and sometimes r' band GROND photometric data (red points) at 
frequencies larger than log(«//Hz) = 15.4 (i.e. after the Lya line). 
The Lya forest is caused by intervening Lya clouds, present in ran- 
dom density and distance along each line of sight. Since the effects 
of this absorption are unlikely to be predicted, we do not try to fit 
the data at frequencies larger than the Lya. [blue and cyan lines in 
Fig. l5UT0lbecome dashed at v larger than log(i^/Hz) = 15.4]. 

Another recurrent feature is the increase in luminosity in the 
redder WISE bands with respect to the decreasing tendency that 
the other IR data have toward lower frequencies. Where the lowest 
frequency WISE data is not an upper limit, indeed, it has always a 
luminosity at least half a decade larger than one would expect ex- 
trapolating the flux from larger frequencies. This could be ascribed 
to the presence of a different component, such as the synchrotron 
hump or an emission from the dusty torus (Calderone, Sbarrato & 
Ghisellini2012). 

The most useful feature that characterizes the optical-NIR 
spectra of our sources is the possibility to directly observe the vL v 
peak of the accretion disc spectrum. Just before the Lya forest ab- 
sorption the SDSS data highlight a peak in the vL v emitted spec- 
trum, that ends the clear increasing trend of the flux from GROND. 
The highest redshift objects show this peak even in the GROND 
bands. This was already pointed out in Sbarrato et al. 2012b, where 
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45.03 


9.17 


45.97 


145212.86 +023526.3 


4.889 


365 


-16.21 
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151002.92 +570243.3 


4.309 


13463 


-14.45 


-12.27 


1.51 


Yo09 


47.08 


44.87 


8.51 


45.82 


152028.14+183556.1 


4.123 


104 


-16.01 








47.07 


44.61 


9.28 


45.43 


165913.23 +210115.8 


4.784 


637 


-15.39 








47.10 


44.12 




45.06 


172026.68 +310431.6 


4.669 


306 


-15.83 








46.96 


44.80 


9.73 


45.74 


172007.19+602824.0 


4.424 


124 


-16.15 








46.98 


44.84 


8.85 


45.79 


222032.50 +002537.5 


4.205 


4521 


-14.87 








46.93 


45.05 


9.15 


46.00 



Table 1. Sources from the DR7 Quasar Catalog with z > 4 and radio-loudness R > 100. Col. [1]: right ascension and declination (i.e. SDSS name: 
"SDSS J. . . "); Col. [2]: redshift; Col. [3]: radio-loudness defined as in SI 1 (R = F$qh z / F j^, fluxes rest frame); Col. [4]: logarithm of the radio flux in 
erg cm -2 s -1 ; Col. [5]: logarithm of the X-ray flux in the energy range 0.3-10keV, in erg cm -2 s — 1 ; Col. [6]: X-ray photon index; Col. [7]: references for 
X-ray data. Ba04 stands for Bassett et al. (2004), Sbl2 for Sbarrato et al. (2012b), Sh05 for Shemmer et al. (2005), Yo09 for Young, Elvis & Risaliti (2009). 
Col. [8]: logarithm of the bolometric luminosity, as obtained by SI 1 ( erg s _1 ); Col. [9]: logarithm of the CIV luminosity, as obtained by SI 1 ( erg s _1 ); Col. 
[10]: logarithm of the BH mass, as derived by SI 1 (solar masses); Col. [1 1]: logarithm of the Broad Line Region luminosity, calculated from the CIV line 
luminosity (ergs -1 ). 



GROND observations allowed to clearly draw the peak of the ther- 
mal emission of B2 1023+25 (z = 5.3). 

The visibility of the accretion disc, along with its vL v peak, 
allow us to fit it with a theoretical model, in order to derive some 
basic properties, such as the overall disc luminosity Ld and the 
black hole mass Mbh- 

4.1 The disc Luminosity 

One way to derive the overall disc luminosity is through the total 
luminosity of the Broad Line Region (BLR). The accretion disc is 
the principal source of UV photons that ionize the BLR, hence the 
two luminosities have to be strictly related. The BLR re-emits a 
fraction C of the ionizing luminosity, that corresponds to the cov- 
ering factor of the BLR with respect to the disc. The average value 
of C is 0.1 (Baldwin & Netzer, 1978; Smith et al. 1981), with a 
large dispersion depending on the geometrical features of the BLR 
itself. The luminosity emitted from the whole BLR can be derived 
from the luminosities of the principal broad emission lines through 
the templates calculated by Francis et al. (1991) and Vanden Berk 



et al. (2001). In this way we obtain a range of possible values of the 
disc luminosity, determined by the uncertainties on the covering 
factor and on the templates. 

On the other hand, when the peak of the accretion disc com- 
ponent is directly visible, we can strongly reduce the uncertainties 
on its bolometric luminosity. This often occurs for the sources in 
our sample, despite the possible contribution of the jet synchrotron 
emission. This is due to the large redshift of the source, red-shifting 
the peak of the disc component within the optical band, red-ward 
of the Lyman-a limit, and because our sources, being very power- 
ful, all have a synchrotron component peaking in the far IR or mm 
band. Hence, with the good coverage of the optical-NIR band given 
by the SDSS, GROND and WISE data, we obtain for all our blazar 
candidates the peak luminosities and can constrain the total disc lu- 
minosity. The uncertainty of our method is hence associated only 
to the measurement errors of the flux and to its possible variability, 
since the GROND, WISE and SDSS data are not simultaneous. 

The upper panel of Fig.[3]shows the direct comparison of our 
estimates of the disc luminosity with half the bolometric luminos- 
ity as defined in Richards et al. (2006). Note that no peculiar trend 



© 0000 RAS, MNRAS 000, 000-000 



Blazar candidates at z > 4 5 



Name 


Date 


Mid time 


Exp: opt/IR 


Seeing 


Mean airmass 




yyyy-mm-dd 


[UTC] 


[s] 


[arcsec] 




QFlQQ 1071047 \fx OOIRIR A 
ctUoct JUZlU47.10-UUlolo.4 


701 7 OS 1 o 


1 o- 1 n-^R 

1U. 1U. jo 


Q1 Q/QAO 

y iy/you 


1 7 
1.7 


i i ^ 

1.10 


QFlQQ T070A77 7 1 j-OOAfi^'} A 
olJSo JU7U47 /.Z1+UU40.JO.0 


701 7 OS 7 1 
ZUiZ-Uo-Zl 


HQ- 1 7-7Q 

uy. iz.zy 


£07/770 

ouz/ / ZU 


o o 
u.y 


1 1 7 

1.1 / 


cnc c toq<i 1 1 An 1 1 A'-ii'xn c 
^Uac> JUojI 1 1.0U+14Z33 /.o 


7m 7 n/i 7n 
ZU 1 Z-U4-ZU 


77 -77-7Q 


yiy/you 


1 7 
1.3 


I.Jo 


cnCC TOQ1 C7/1 1C i n£1&<1 7 

:>D:>c> JUy loZ4.3o+Uo3o:>3. J 


7H1 7 D/1 7 1 
ZU 1 Z-U4-Z 1 


77 -77 -nn 
Ld.LI .UU 


yiy/you 


1 A 
1.4 


1 7C 

I.Zj 


CT^CC T1AO*tT2 A7_i_7^/17<Q A 

oDIjo J lUZoZ3.oZ+Z74Z:>y.O 


701 7 CiA 1 A 
ZU I Z-U4- 1 


7 7 ■ 7 O ■ 1 7 

Z7.zy.i7 


yiy/you 


1 7 
l.Z 


7 1 1 
Z. 1 1 


CT^C C T1 HQ 7 1 7 77i1Q77H7 1 

oDoo J 1UJ / 1 /. / J+18Z3U3. 1 


701 7 OA H/1 

ZU 1Z-UO-U4 


77-1 Q ■ 7 7 


yiy/you 


1 A 
1.4 


1.D4 


oDijo J !Uj3ZU.4Z-UU104y. j 


701 7 (\£L HQ 

ZUlZ-Uo-Uo 


77 -77-^^ 
Lj.LL. J D 


yiy/you 


1 n 
l.U 


1 io 
1.1(5 


^Uac> J 1Z3jU3.UZ-UUU331.0 


701 7 HA 1 /I 

ZUlZ-Uo-14 


7 7 ■ 1 n ■ /I n 
ZJ.1U.4U 


oU3/o4U 


1 7 

1.3 


1.10 


MJdS J13U/3o.o3+lMJ/:>z.l 


ZU 1Z-U/-U4 


77. 1 n 

23:1 /:3/ 


yiy/you 


7 A 

Z.o 


1 /i 1 
1.41 


cnCC TIQHII 77_i_77777C A 
^LJa^ J 13 1 1Z1.3Z+ZZZ /3o.O 


7H1 7 OA 7H 

ZU1Z-U0-3U 


Ul .U/.4Z 


yiy/you 


1 1 
1.1 


1 QR 
l.OO 


QFlQQ T1 77^ 1 7 AQ_i_1 1 7770 7 

oUoo j i jzj iz.4y+i izjzy. / 


701 7 07 OA 
ZU1Z-U /-U0 


77-OA-77 
Z7.U0.Z7 


Q70/QAO 

yzu/you 


1 1 
1 . 1 


1 77 
1.3Z 


SDSS J141209.96+062406.8 


2012-08-26 


23:27:05 
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1.75 


SDSS J142048.01+120546.0 


2012-08-12 


23:34:17 


919/960 


2.9 


1.62 


SDSS J143003.96+144354.8 


2012-08-13 


00:13:27 
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2.0 


1.82 


SDSS J143413.05+162852.7 


2012-08-21 


23:31:37 


672/840 


1.8 


1.90 


SDSS J145212.85+023526.4 


2012-08-28 


23:33:18 


919/960 


1.5 


1.23 


SDSS J152028. 14+183556.1 


2012-09-04 


23:38:49 


919/960 


1.8 


2.04 


SDSS J165913. 23+210115. 8 


2012-09-15 
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2012-05-28 
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Table 3. Log of the GROND observations. Exposures refer to optical/NIR filters while the average seeing is calculated in the r' band. 
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19 


,17±0.02 


19 


,28±0 


05 


19 


,13±0 


.07 


18 


.96±0. 


14 


0.080 


SDSS J030437 


21+004653.6 
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,34±0.04 
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,32±0.04 


20 
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SDSS J085111 
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10 
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SDSS J105320. 
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03 


19 
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.35±0.01 
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31 


21 
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15 
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20 
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69 


0.062 
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08 


19 
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20 
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,56±0. 
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Table 4. Observed magnitudes, not corrected for Galactic foreground reddening, in the AB system. Errors include systematics. The last column reports the 
value of the Galactic Av from Schlafly & Finkbeiner (2011). 



nor an evident offset is present between the two luminosities, but 
they are clearly linearly correlated. This represents a further confir- 
mation of the "naked" accretion disc. Indeed, the bolometric lumi- 
nosity in Richards et al. (2006) is derived from the observed flux 
measured at specific observed wavelengths. If the flux was a com- 
bination of the emissions from the accretion disc and an additional 
component, e.g. the synchrotron emission, it would not be a good 
tracer of the accretion disc emission. In that case, it would not be 
so tightly correlated with our estimate of the overall accretion disc 
luminosity, that suffers less from possible contaminations, being 
directly derived from the peak of the disc emission. 



4.2 Black Hole Mass Estimate 

The black hole mass is a key feature to describe the quasars in our 
sample. The high redshift, however, introduces some difficulties in 
applying the commonly used virial method. Indeed, the CIV line is 
present only for z < 4.7 in the SDSS spectral range, hence it is not 
possible to obtain the virial masses for the whole sample. Moreover, 
the CIV line is the less reliable broad emission line to apply the 
virial argument (Shen 2013 and references therein). Therefore, to 
derive the Mbh coherently for all our objects, we apply a different 
method. 
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As mentioned above, the extreme power of these objects en- 
sure that the synchrotron hump is at frequencies low enough to 
leave "naked" the spectrum produced by the accretion disc, that 
peaks in the optical-IR band due to the very high redshift. Hence, 
with a good coverage of this wavelength range, the accretion disc 
emission can be directly modeled and fitted. Since bright broad Lya 
and CIV lines are visible in the optical spectra of our sample, we 
can reasonably assume that the accretion disc is radiatively effi- 
cient, being able to ionize the broad line region (BLR). We thus 
apply the simplest hypothesis for a radiatively efficient accretion 
disc, i.e. the Shakura-Sunyaev disc (Shakura & Sunyaev, 1973). 
According to this model, the geometrically thin, optically thick ac- 
cretion disc can be fitted with a multicolor black body spectrum. 
This spectrum depends only on the central black hole mass (A/bh) 
and on the mass accretion rate (A/), that is directly traced by the to- 
tal luminosity emitted by the disc (Ld = r\M<? , with rj = O.OcQ). 
In our case, we can precisely constrain both parameters. 

The same model has been applied by Calderone et al. (2013) 
to fit the accretion discs of a sample of radio-loud Narrow Line 
Seyfert Is, and in that paper the model is described in detail. As we 
already pointed out, at these redshifts, we clearly see the vL v peak 
of the accretion disc component. This allows to directly derive the 
disc luminosity Ld of all our sources [i.e. Ld = 2v p L Vp , where v p 
is the peak frequency of the disc spectrum]. 

With the total disc luminosity fixed by the vL v peak, we are 
left with only the black hole mass as a free model parameter. In 
a multicolor blackbody model, for a fixed value of the disc lumi- 
nosity, the black hole mass affects only the peak frequency, since it 
determines the emitting surface and its temperature. In other words, 
a larger black hole mass implies a larger disc surface, that needs to 
be colder to emit a fixed Ld. This implies a rigid shift of the vL v 
peak towards smaller frequencies (see also Fig. 3 in Calderone et 
al. 2013). Therefore, if the overall disc luminosity is determined 
through the i/ p L„ p , it is straightforward to estimate the central 
black hole mass, with a good sampling of the optical-NIR energy 
range. 

We have derived the black hole mass for all the sources in 
our sample. Figs.lSlflOlshow the SEDs of our objects together with 
three models with same disc luminosity but different masses. Note 
that the models do not fit the data at log v > 15.4, i.e. in the Lya 
forest region, since the effect of the absorption produced by ran- 
domly distributed Lya clouds along different lines of sight is un- 
likely to be predicted. The lowest and highest masses indicate the 
uncertainty associated to the measurement errors, possibe presence 
of other contributions, and non strict simultaneity of the data. The 
best values of the black hole masses we derived are listed in Tab. 
[2] The indicated errors are not formal errors, but the mass range 
determined by applying the three models to the same data. 

Since our sample is composed of highly luminous sources, 
given the survey flux limit and the high redshift, we expect a large 
black hole mass, if the objects emit at some large Eddington frac- 
tion, but do not exceed the Eddington limit. This is indeed what 
we find: all quasars in our sample host black holes with masses 
Mbh > 10 9 Mq. The red histogram in Fig. [2] shows the distri- 
bution of masses that we obtain for our sample using the fitting 
method. Our results are distributed on a lognormal distribution, 



The efficiency can be defined as rj = GM B - tl /(2R iu c 2 ), where R in 
is the inner radius of the accretion disc, for which we assume a value of 
ili n = 6GA/bh/ c2 . as in the case of a non-rotating black hole. With this 
simple hypothesis, we assume r\ = 0.08. 
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Figure 2. Normalized distribution of the black hole masses of the quasars 
in our sample. In blue, the histogram of the results obtained by SI 1 in their 
work. In red, our results. The dashed lines are best fit log-normal distribu- 
tions, with values: (log M^h/Mq) = 9.37; a = 0.43 (blue histogram, 
results of Sll), (log M-bh/Mq) = 9.31; a = 0.21 (red histogram, this 
work). 



with an average value: 

M BH 



log 



Mr. 



= 9.31 ±0.21 



(2) 



The blue histogram shown in Fig. ^corresponds to the values of the 
black hole mass derived by S 1 1 through the virial method (namely, 
through the FWHM of the CIV line). 



los 



Mbh \ 

Mq / Sll 



9.37 ±0.43 



(3) 



Comparing the two distributions, despite the paucity of objects, we 
note that the average value is the same, but the dispersion obtained 
with the accretion disc fitting method is smaller. The difference in 
dispersion is evident from the lower panel of Fig. [3] too. This plot 
compares directly our Mbh estimates with the virial results (we in- 
clude only the 14 sources for which Sll were able to virially derive 
the A/bh). The large dispersion of the Sll results prevents us from 
a clear "source by source" comparison, but it is reasonable to assess 
that no systematic error seems to affect our estimates with respect 
to the virial black hole mass estimates. Nevertheless, the great dif- 
ference between the dispersions of the two methods is clearly an 
interesting issue, that will be investigated more deeply in the next 
section. 



5 DISCUSSION 

In this work we constructed a sample of extremely radio-loud 
high-redshift quasars, in order to select a group of good blazar 
candidates. We studied their thermal emission, deriving the lu- 
minosity emitted by the accretion disc (assumed as an optically 
thick, geometrically thin Shakura-Sunyaev disc) and the central 
black hole mass. We took advantage of the good data coverage that 
can be achieved combining the WISE, SDSS and GROND data. 
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Figure 3. Compaiison of Sll with our results, including only the 14 
sources out of our sample of 19 candidates with both the virial black hole 
mass and the accretion disc fitting estimates. Upper panel: the disc lu- 
minosity we obtain with the accretion disc fitting compared to half the 
bolometric luminosity (Calderone et al. 2013) as derived by Sll. SDSS 
J102623.61+254259.5 (B2 1023+25) is not included in this plot, since its 
redshift excludes from the SDSS wavelength window all the wavelengths 
from which the bolometric luminosity is typically derived (Richards et al. 
2006). Lower panel: Values of the black hole masses as derived in this work 
with the accretion disc fitting (along with our error), compared with the 
virial masses derived by Sll. This plot shows only the 14 objects with red- 
shift that allows the virial black hole mass estimate. The 5 candidates with 
the highest redshift, do not show the CIV broad line int the SDSS wave- 
length interval. We assigned to the virial masses an average uncertainty of a 
factor 3.5 (Vestergaard & Peterson, 2006). Note that the distribution of our 
.\/. ; . I estimates is definitely less dispersed than the distribution of the virial 
results (as can be clearly seen in Fig. [2] Our results show no systematics 
with respect to the Sll results. 



The optical-IR wavelength range, indeed, is where the Shakura- 
Sunyaev spectrum is expected to peak for SMBHs at z > 4. The 
peak of the emitted spectrum is indeed visible for our objects, and 
this greatly help to obtain good fits with a Shakura-Sunyaev disc 
model. 

It is worth to notice that, once we assume an accretion disc 
model, the uncertainty of this method is strictly linked to the data 
uncertainties (i.e. the larger the data errors, the less constrained our 
estimate). The wavelength coverage that we have in the IR-optical- 
UV band and the good precision of the available data, allow us to 
achieve an uncertainty of a factor of two. Note that, in principle, 
we can obtain even more precise results with more accurate data. 
In the virial method, instead, the precision of the mass estimates 
are driven by the uncertainties of the calibrating model, that can- 
not be reduced below a factor 3-4 (Vestergaard & Peterson, 2006). 
Furthermore, the FWHM of the broad emission lines could be af- 
fected by two systematic effects. The first is related to the geometry 



of the BLR. If the clouds forming the BLR are distributed in a flat- 
tened geometry, their FWHM would be viewing angle dependent, 
being smaller for small viewing angles (Decarli, Dotti & Treves, 
201 1 ). In turn, small viewing angles are expected for sources with 
a large radio-loudness. Therefore, if such an effect plays a role, 
our sources would have virial black hole mass estimates that are 
systematically smaller than what derived from the same sources 
viewed at larger angles. On the other hand, our accretion disc fit- 
ting is much less viewing angle dependent, and thus the comparison 
between the two black hole mass distributions can tell something 
about the BLR geometry. 

The second effect concerns the effect of radiation pressure 
upon the line emitting clouds, as discussed by Marconi et al. (2008, 
2009). Accounting for this extra force makes the virial estimate of 
the black hole mass larger, and more so for near-Eddington disc 
luminosities. Again, if the Shakura-Sunyaev disc is a reasonable 
description of the reality even at luminosities near Eddington, then 
the estimate of the black hole mass is independent of the disc to 
Eddington luminosity ratio Ld/L^dd- 

The comparison of the two black hole mass distributions gives 
these important results: 

(i) The average black hole mass is the same — As explained 
above, this suggests that the effect of both a flattened distribution 
of the clouds and/or the radiation pressure of the disc radiation play 
a small role. What inhibits a stronger conclusion is only the small 
number of objects, but this issue surely deserves some attention in 
the future, when larger samples will be available. 

(ii) The intrinsic width of the distribution must be extremely nar- 
row — As explained above, the virial method relies on calibrating 
some scaling relation, as the size-luminosity relation, whose intrin- 
sic scatter does not allow a determination of the black hole mass 
more accurate than a factor 3 or 4. The disc fitting method, instead, 
does not suffer from these limitations, but depends only on the pre- 
cision of the data. For our sources, on average, we can assign an 
uncertainty of a factor 2. Fig. [2] shows indeed that our results are 
less dispersed than the results obtained by S 1 1 . 

Each observed dispersion (cr bs) is actually the convolution of 
the intrinsic dispersion of the black hole masses (o\m bh ) w ' m me 
typical error that is done during each of the two measurements 

(o"err): 



2 

°obs 



(4) 



Therefore, with two different observed distributions, we could be 
able to estimate the intrinsic dispersion of masses in our sample. 

Bearing in mind that the sample is composed by only 19 ob- 
jects, and that Sll obtained virial black hole masses for only 14 
of these sources, let us compare the observed dispersions obtained 
by Sll (a^l = 0.43) with our dispersion (cr^; work = 0.21). 
They are of the same order of the typical errors of the two methods 
(a^ 1 ~ 0.4 - 0.5 and cr^ s wo,k ~ 0.2 - 0.3). This suggests an ex- 
tremely narrow intrinsic distribution, almost a 5-function. Interest- 
ingly, a very narrow Mbh distribution is also found by Calderone 
et al. (2013) for a sample of NLS1 galaxies, even if they consider a 
different range of masses and features of the sample. Nevertheless, 
this result has to be considered carefully. First of all, we still have 
too few objects to draw statistically relevant conclusions. More- 
over, it could be affected by selection effects. Indeed, our selec- 
tion criteria collect the most luminous radio-loud quasars from the 
SDSS, that may already be composed by the "tip of the iceberg" 
of the quasar population at high redshift. Since it is reasonable to 
expect that the most luminous are also the most massive quasars, 
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Figure 4. Disc luminosity as a function of the black hole mass, both in 
logarithmic scale. Blue squares are the results of Sll for the objects in our 
sample. Specifically, the Mbh is obtained with virial methods, while the 
is half the bolometric luminosity derived as in Richards et al. (2006) 
(Calderone et al. 2013). Red dots are the results of the method used in this 
work for the 19 objects in our sample. Yellow dots are all the radio-quiet 
AGN included in Sll, with L d = L bo i/2. The grey stripe highlights the 
minimum luminosity measured by S 1 1 among the radio-quiet AGN. The 
dashed lines indicate accretion regimes at the Eddington rate and at a hun- 
dredth of Eddington rate, as labelled. Note that Sll obtain results that goes 
at disc luminosities larger than the Eddington limit. 



our sample likely includes the objects with the most extreme black 
hole masses. This possible selection effect may partially justify the 
narrowness of the intrinsic distribution, if we preferentially select 
the most luminous and massive objects. Bearing in mind these ob- 
servations, we can conclude that our small sample leads to the sug- 
gestion that the extremely radio-loud quasars located at z > 4, and 
present in the SDSS, tend to have the same black hole mass, always 
larger than 10 9 M Q . 

We now turn our attention to the comparison between the disc 
properties of radio-loud and radio-quiet AGN at z > 4: 

(i) The thermal emission of radio-loud AGN do not differ from 
the analogous thermal emission of radio-quiet AGN — Fig. [4] 
shows the overall disc luminosity as a function of the black hole 
mass for the objects in our sample, as measured both by SI 1 (blue 
squares) and with our method (red dots), compared with all the 
radio-quiet AGN studied by Sll (yellow dots) above z — 4. For 
the data obtained by Sll, to estimate the disc luminosity we used 
half the bolometric luminosity calculated in their work, following 
Richards et al. (2006) and Calderone et al. (2013). Fig.|4]shows that 
there are no intrinsic differences between radio-loud and radio- 
quiet AGN in their black hole mass and accretion luminosities. 
Therefore the thermal emission from the accretion disc and the 
mass of the central black hole of AGN, in presence or absence of a 
jet, do not differ. 

(ii) Extremely radio— loud AGN at high redshift have high accre- 
tion rates — Fig.[4]illustrates also the range of accretion rates of the 
objects in our sample, in units of the Eddington one. The two grey 
dashed lines in fact correspond to a Eddington rate (L c \ ~ -I/Edd) 



and a hundredth of Eddington (id = 10~ I/Edd). They are gener- 
ally considered as the two physical limits of a radiatively efficient 
accretion operated by a Shakura-Sunyaev accretion disc. Note that 
two extremely radio-loud objects, with the data obtained by Sll, 
have a super-Eddington accretion rate. On the contrary, our results 
never exceed this limit. This suggests that at least in these cases 
our method is more reliable. It is also clear from Fig. |4]that all our 
sources have high accretion rates. Thus our sources, besides having 
a relatively large black hole mass, [log(JWBH/A^e) = 9.3], accrete 
at ~10% Eddington, on average. High rates of accretion justify the 
presence of many extremely massive black holes at high redshift, 
that perhaps are just ending the "fast growing age" through accre- 
tion at the Eddington rate. Indeed, the formation and growth of the 
first supermassive black hole is an open issue, that our results will 
hopefully help to solve. 

(iii) Blazar candidates — The sample, selected through a cut in 
radio-loudness, preferentially collects jets aligned with the line of 
sight. Among them, we expect to find a number of blazars, namely 
sources with a viewing angle comparable to or smaller than the 
beaming angle 1 /Y. If we knew the viewing angle of our sources, 
we could reconstruct the number density of the entire (parent) pop- 
ulation. Specifically, if a source has a viewing angle smaller than 
the beaming angle, with an estimate of Y (through modelling and/or 
superluminal motion) we would be able to derive the total num- 
ber of analogous objects, This is indeed the final goal of our re- 
search. We have already identified as a blazar the more distant of 
the objects in our sample: B2 1023+25 at z ~ 5.3 (Sbarrato et 
al. 2012b). In our previous work, we suggested a viewing angle 
smaller than 1 /Y for this source, as it is for the other two blazars 
of the original sample (SDSS J083946.22+5 11202.8 and SDSS 
J151002.92+570243.3, see jS), Each of these blazars implies the 
existence of other ~ 450(r/15) 2 misaligned sources sharing the 
same characteristics, including the presence of an actively accret- 
ing supermassive black hole, in the same region of the sky covered 
by the SDSS+FIRST survey. 



6 CONCLUSIONS 

In this work, we composed a sample of 31 extremely radio-loud 
quasars at z > 4. We estimated the black hole masses of the 19 
objects among them that are visible from La Silla (Chile), and that 
are hence observable with GROND. The GROND data allowed us 
to observe directly the peak of the accretion disc spectrum. Hence, 
combining them with WISE and SDSS data, we have been able to 
study the thermal emission of those objects, and to derive with a 
precision of a factor of two their black hole masses. We found a 
very small range of masses for our sample, that peaks at Mbh = 
1O 93 M . 

Our method, along with the good data coverage of the optical- 
IR band, allowed us to derive both the black hole mass and the 
accretion disc luminosity of the objects in our sample. Hence we 
were able to estimate the accretion rate of the 19 extremely radio- 
loud quasars. We found high accretion rate values for all of them. 
This could introduce interesting hints on the formation and (possi- 
bly fast) growth of the first supermassive black holes at very high 
redshift. 

Our sample collects highly radio-loud quasars with extremely 
massive black holes. Such a radio-loudness criterion preferentially 
selects jets aligned near our line of sight. We expect to find in 
our sample a number of blazars, i.e. quasars seen with a viewing 
angle comparable to (or smaller than) the jet beaming angle. We 
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will derive the jet orientation of the best candidates by fitting their 
broad-band SEDs. Specifically, we have planned for our three best 
candidates a set of X-ray observations with the Swift X-ray Tele- 
scope. X-ray data will provide hints on the jet non-thermal emis- 
sion, helping us to constrain the jet viewing angle. 

Thanks to their peculiar orientation, by identifying these ob- 
jects as blazars we would be able to derive the total number of 
analogous extremely massive radio-loud objects with jet directed 
in other directions. In other words, with our approach we expect 
to obtain the overall number of the parent population of extremely 
massive blazars in the early Universe. This would allow a direct 
comparison of the density functions of highly massive radio-quiet 
and radio-loud objects even at very high redshift. 
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Figure 5. IR-optical-UV SEDs of the quasars observed by GROND in 
the rest frame, together with standard accretion models. In each figure, the 
black line is the SDSS spectrum. Red and magenta points are data from 
GROND and WISE, respectively. Green points are archival data (ASDC 
SED builder, NED). The grey stripe indicates the vL v peak luminosity of 
the disc. The blue line is the multicolor black body that best fits our data, 
calculated with L d and A/bh values reported in each figure. Outside the 
range marked by the cyan lines, the disc model cannot fit satisfactory the 
data. Note that at log v > 15.4 the models do not fit the data, because of 
the Lya forest. 
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Figure 6. See caption of Fig. [5] 



Figure 7. See caption of Fig. \5\ 
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Figure 8. See caption of Fig. [5] 



Figure 9. See caption of Fig. \5\ 
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Figure 10. See caption of Fig.[5] 
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